home *** CD-ROM | disk | FTP | other *** search
/ Acorn User: China / Acorn User China CD-ROM (UK) (Disc B) / Acorn User China CD-ROM (UK) (Disc B).bin / STUTTGART / LANG / GOFER / scripts_GaussInt < prev    next >
Encoding:
Text File  |  1989-03-31  |  1.7 KB  |  47 lines
  1. {------------------------------------------------------------------------
  2.             Factorizing Gaussian integers using Gofer               
  3. ------------------------------------------------------------------------}
  4.  
  5. prime_divisors :: Int -> [Int]
  6. prime_divisors x = f (abs x) where
  7.   f (2*n) = 2:f n                             -- deal with even case
  8.   f (3*n) = 3:f n
  9.   f (5*n) = 5:f n
  10.   f (7*n) = 7:f n
  11.   f n | n<2       = []
  12.       | otherwise = d:f (n/d) where
  13.             d                         = maybe 11
  14.             maybe k | n `mod` k == 0  = k
  15.                     | k*k > n         = n            -- n has to be prime
  16.                     | otherwise       = maybe (k+2)  -- next odd candidate
  17.  
  18.  
  19. sum_of_squares :: Int -> [Gauss Int]
  20. sum_of_squares n = [ (a:+!b) | a<-[1..m], b<-[1..(n-a*a)], n == a*a+b*b]
  21.             where m:_ = [ k-1 | k <- [1..], k*k > n ]
  22.  
  23. irreducible :: Int -> [Gauss Int]
  24. irreducible p | p == 2         = [unit + i, unit - i]
  25.               | p `rem` 4 == 3 = let x = p:+!0 in [x,x] 
  26.               | otherwise      = sum_of_squares p
  27.  
  28. divides :: (Gauss Int) -> (Gauss Int) -> Bool
  29. z `divides` z' = a `rem` d == 0 && b `rem` d == 0 where  
  30.                 a :+! b  = z' * (conjugate z)                   
  31.                 d       = norm z
  32.  
  33. factors :: (Gauss Int) -> [Gauss Int]
  34. factors z | z == zero  = []
  35.           | norm z == 1  = [z]
  36.           | otherwise    = w:factors (z/w) where
  37.               w | u `divides` z = u 
  38.                 | otherwise     = u'
  39.                            where
  40.                             u:u':_ = irreducible p
  41.                             p:_    = prime_divisors (norm z) 
  42.  
  43. factorise :: Gauss Int -> String
  44. factorise = show . factors 
  45.  
  46. ----- End Complex
  47.